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ABSTRACT 

We discuss three modes of oscillation of accretion disks around rotating magnetized neutron stars 
which may explain the separations of the kilo-Hertz -periodic oscillations (QPO) seen in low mass X- 
ray binaries. The existence of these compressible, non-barotropic magnetohydrodynamic (MHD) modes 
requires that there be a maximum in the angular velocity fi^(r) of the accreting material larger than 
the angular velocity of the star f2*, and that the fluid is in approximately circular motion near this 
maximum rather than moving rapidly towards the star or out of the disk plane into funnel flows. Our 
MHD simulations show this type of flow and Q,t,(r) profile. The first mode is a Rossby wave instability 
(RWI) mode which is radially trapped in the vicinity of the maximum of a key function g(r)!F{r) at tr. 
The real part of the angular frequency of the mode is u> r = mf^(rft), where m = 1, 2... is the azimuthal 
mode number. The second mode, is a mode driven by the rotating, non-axisymmetric component of 
the star's magnetic field. It has an angular frequency equal to the star's angular rotation rate f2*. This 
mode is strongly excited near the radius of the Lindblad resonance which is slightly outside of 7\r. The 
third mode arises naturally from the interaction of flow perturbation with the rotating non-axisymmetric 
component of the star's magnetic field. It has an angular frequency f2»/2. We suggest that the first 
mode with m = 1 is associated with the upper QPO frequency, v u ; that the nonlinear interaction of the 
first and second modes gives the lower QPO frequency, vt = v u — v*] and that the nonlinear interaction 
of the first and third modes gives the lower QPO frequency vi = v u — v„/2, where 2/* = 0,^/2-k. 

Subject headings: keywords: accretion, accretion disks — stars: neutron — X-rays: binaries — 
magnetohydrodynamics — Instabilities — Waves 



1. INTRODUCTION 

Low mass X-ray binaries often display twin kilo-Hertz 
quasi-periodic oscillations (QPOs) in their X-ray emissions 
(van der Klis 2006; Zhang et al. 2006). A wide variety of 
different models have been proposed to explain the ori- 
gin and correlations of the different QPOs. These include 
the beat frequency model (Miller, Lamb, & Psaltis 1998; 
Lamb & Miller 2001; Lamb & Miller 2003), the relativistic 
precession model (Stella & Vietri 1999), the Alfven wave 
model (Zhang 2004), and warped disk models (Shirakawa 
& Lai 2002; Kato 2004). 

A puzzling aspect of the some of the twin QPO sources 
considered in this work is that the difference between 
the upper v u and lower vi QPO frequencies is roughly 
cither the spin frequency of the star v* (3 cases where 

= 270, 330, & 363 Hz) or one-half this frequency, v*/2 
(4 cases where v* = 401, 524, 581, & 619 Hz), for the 
cases where is known, even though v u and vi vary sig- 
nificantly (see, e.g., Zhang et al. 2006). A further type 
of behavior appears in the source Cir X-l (Boutloukos et 
al. 2006), but this is not considered here. The cases where 
Vu~ v i ~ v * may be explained by the beat frequency model 
(Miller et al. 1998), but the explanation of the cases where 
v u — vu ~ f*/2 is obscure. 

Li and Narayan (2004) analyzed the stability of an in- 
compressible rotating flow where the radial profiles of the 
plasma angular rotation rate O^, the magnetic field B z , 
and the density p change sharply at the magnetopause ra- 
dius but are independent of z. The magnetic field was 
dynamically important in the sense that pu 2 ~ B 2 /47r 
with u = ril^cb . They found Rayleigh- Taylor and Kelvin- 
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Helmholtz type instabilities in the vicinity of the magne- 
topause. Recently, Fu and Lai (2009) have done a system- 
atic analysis of short wavelength modes of a disk for cases 
where the magnetic field is not dynamically important. 

The present work is a continuation of the study of 
Lovelace and Romanova (2007; hereafter LR07) of the 
magnetohydrodynamic (MHD) stability of the compress- 
ible, non-barotropic stability of the boundary between an 
accretion disk and the magnetosphere of a rotating mag- 
netized star for conditions where S7<^, B Zl and p vary 
smoothly with radial distance but are independent of z, 
and where the magnetic field is dynamically important. 
The radial profiles of these quantities are known for dif- 
ferent conditions from MHD simulations (Romanova et 
al. 2002; Long, Romanova, & Lovelace 2005; Romanova, 
Kulkarni, & Lovelace 2008; Kulkarni & Romanova 2008). 
The recent MHD studies (Romanova et al. 2008; Kulkarni 
& Romanova 2008) discovered conditions where a global 
Rayleigh- Tayor instability occurs and changes the nature 
of the accretion flow from a regular funnel flow pattern 
to a chaotic flow of plasma fingers. The present work is 
concerned with a radially localized instablity of a globally 
stable approximately axisymmetric flow. Of particular im- 
portance to the present work is that tt^r) is observed to go 
through a maximum significantly larger than the angular 
rotation rate of the star $1* = 2~kv*. The importance of the 
maximum of O^(r) for models of QPOs was discussed ear- 
lier by Alpar and Psaltis (2005). LR07 showed that there 
was a Rossby wave instability (RWI) or unstable corota- 
tion mode with azimuthal mode number m = 1,2.. radially 
trapped at tr near the maximum of O^(r). It was sug- 
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gested that the upper QPO frequency was ui r = mf^r^) 
and that the lower QPO frequency was due to the inter- 
action of the mode with rotating non-axisymmetric field 
of the star. The theory of the RWI was developed by 
Lovelace et al. (1999) and Li et al. (2000) for accretion 
disks and earlier by Lovelace & Hohlfeld (1978) for disk 
galaxies. Hydrodynamic simulations of the RWI instabil- 
ity in disks were done by Li et al. (2001), while Sellwood 
& Kahn (1991) used N— body simulations to study the 
instability in galaxies. The instability has important role 
in the accretion-ejection instability of disks discussed by 
Tagger and collaborators (e.g., Tagger & Varniere 2006; 
Tagger & Pellat 1999). 

Section 2.1 develops the general compressible, non- 
barotropic MHD equations for free perturbations an ax- 
isymmetric equilibrium flow with no z— dependence. Sec- 
tion 2.2 treats the driven perturbations due to the ro- 
tating non-axisymmetric component of the star's mag- 
netic field. Section 2.3 treats the driven-modulated per- 
turbations which result from the interaction of the non- 
axisymmetric field of the star and the flow perturbation. 
Section 3 develops a detailed model of the axisymmetric 
flow/field equilibrium based on results from MHD simu- 
lations. Section 4 discusses the results obtained applying 
the theory of §2 to the model of §3. Section 4.1 treats the 
free perturbations, §4.2 the driven perturbations, and §4.3 
the driven-modulated perturbations. Section 5 briefly dis- 
cusses the nonlinear effect of the unstable mode. Section 
6 gives the conclusions of this work. 



2. THEORY 

We consider the stability of the magnetopause of a ro- 
tating star with an aligned dipole magnetic field. The 
envisioned geometry is shown in Figure 1. We use an 
inertial cylindrical (r, 0, z) coordinate system. The equi- 
librium has (d/dt — 0) and (d/dcf) — 0), with the flow 
velocity u = u$(r)<f> = rfl^,(r)(j) . That is, the ac- 
cretion velocity u r and the vertical velocity u z are as- 
sumed negligible compared with u$. The equilibrium mag- 
netic field is B — B(r)z. The equilibrium flow satisfies 
pr(Q 2 K — QV) = —d[p + B 2 /(8Tr)]/dr, where p is the den- 
sity, p the pressure, $ the gravitational potential, and £1% 
the Keplerian angular rotation rate of a single particle. 
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Fig. 1. — Sketch of the envisioned disk, magnetic field, and rotat- 
ing star geometry. Q„ is the angular rotation rate of the star, the 
wiggly lines indicates radiation from the surface of the disk where 
the optical depth is unity, and ff indicates the funnel flow. The re- 
gion outside the disk and funnel flow is occupied by a low-density, 
high-temperature corona. 



2.1. Free Perturbations 

The perturbed quantities are: the density, p = p + 
Sp(r, <f>, t); the pressure is p = p + 5p(r, <fi, t); and the flow 
velocity is u = u+<5u(r, cj), t) with Su = (Su r , Su^, 0). Also, 

B — zB + zSB. The equations for the perturbed flow are 
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where D/Dt = d/dt + u-V and we refer to S = p/ip) 1 as 
the entropy of the disk matter. 

We consider perturbations ~ /(r)exp(im0— iu>t), where 
m = 0, 1, 2, .. is the azimuthal mode number and u> the an- 
gular frequency. For free perturbations u> = uj r + iu>i and 
for the growing modes of interest u>i > 0. (For the driven 
perturbations considered in §2.2 and 2.3, w, = 0.) From 
equation (la), we have 

iAuj Sp = V ■ (p Su) , (2) 

where Aw(r) = lj — mCl$(r) and Q<a = u^/r. 
From equation (lb) we have 
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[r 3 d(r 4 n'j > ) /dr] i is the radial epicyclic fre- 
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From equation (lc) and (Id), we have 
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Combining these expressions, 
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are the length-scales of the entropy S = p/p 1 and B/p 
variations. Also, 

1/2 |_g, 

' = I — I , and <• i = —p== , (56) 



are the sound and Alfven speeds respectively. We denote 
Cf = (c 2 s + c^) 1 / 2 as the fast magnetosonic speed. 
It is useful to introduce 

Sp* 



P 

Equations (3) can then be written as 
ASu r + Bdu^ = f r , 
C 5u<j, + D5u r = , 
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where 
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For a strong magnetic field ca S> c s we have L* — > For 
a weak magnetic field <C c s , we have L» — » Ls- Using 
the equation for the equilibrium pr(Vl 2 K — 0|) = — dp^/dr 
we have 

i r(n 2 K -n 2 ) i 

T" = ~2 - T ' ( 8c ) 



which is useful later. 

From here on we simplify the analysis by neglecting the 
self-gravity of the disk. That is, we neglect the <5<£> terms 
in f r and fa. 

We solve equations (7a) and (7b) to obtain 
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is the 'Lindblad factor.' For real frequencies u>, L)(r) is 
equal to zero at a Lindblad resonance at n,. On the other 
hand at a corotation resonance !K[Aw(r)] = at rp. Using 
equation (8c) we find 
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Equations (9) can now be used to obtain 
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It is useful to note that 
dAtu 



Alu ^ dp^/dr 

_ Mn?-4fi|) 

From equations (4) and (9a) we have 
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Equations (11) and (12) can now be substituted into equa- 
tion (2). When this is done we find that all of the terms 
involving d^ / dr apart from the overbraced term in equa- 
tion (11) cancel out. 

Combining equations (9) with equation (2) gives 
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where g = exp(2 f dr/L*). In the limit of no magnetic field 
(ca — > and L* — > Lg), this equation becomes the same 
as the equation for the Rossby wave instability (Lovelace 
et al. 1999; Li et al. 2000). 

A quadratic form can be obtained by multiplying equa- 
tion (13) by (the complex conjugate of ^) and inte- 
grating over the disk. Assuming r^* (d^ / dr)T / 'Q^ — > 
for r — > 0, oo, we obtain 
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For corotation modes where \Auj\ 2 <C Q|, the third and 
fourth integrals on the right-hand side of equation (14) are 
possibly important. Earlier work by Lovelace and Hohlfcld 
(1978), Lovelace et al. (1999), and Li et al. (2000) found 
that a corotation instability was possible if the quantity 
gT has a maximum or minimum at the radial distance tr 
where 3?(Aw) = 0. For a weak magnetic field {c 2 A <C c 2 
where L* ss Ls), we find gT = S 2 ^T which agrees with 
the result of Lovelace et al. (1999) and Li et al. (2000) 
where the self-gravity of the disk was neglected. With 
negligible self-gravity, instability is found only for condi- 
tions where gT has a maximum as a function of r. For a 
maximum the radial group velocity of the Rossby waves is 
directed towards tr (Lovelace et al. 1999). In the strong 
B-field limit (c 2 A > c 2 where L* w L B ), gT = (B/p) 2 T. 

We can rewrite equation (13) in a form more amenable 
to numerical analysis. That is, 



1 V \ ° Auo + (Aw) 2 
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where V = Vjirp) = (fl^/^T^ 1 , where V is given by 
equation (10b), and T by equation (10a). The primes de- 
note radial derivatives. Also, 
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where g is defined below equation (13). 

If we let * = (V) 1 2 f, then equation (15) can be writ- 
ten in the form of a Schrodinger equation, 
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where U (r) is an effective potential. Both U and ip are in 
general complex for complex ui. 

A quadratic form analogous to equation (14) can be ob- 
tained by multiplying equation (17) by ip* dr and integrat- 
ing over r. Assuming p*p' vanishes at small and large r, 
we obtain 
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The imaginary part of this equation gives 
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where the ellipsis allow for contributions from the Co and 
C2 terms. If these terms are negligible as found in our nu- 
merical evaluations, it follows that a non-zero growth rate 
uji is possible only for conditions where [ln(<jLF)]' changes 
sign as mention above (Lovelace & Hohlfeld 1978). Thus, 
gT is a key function for stability of the considered flow. 

Using equations (4), (9a), and (12) we obtain the rela- 
tion of the temperature perturbation to *, 
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This equation shows that the temperature perturbation is 
strongly enhanced at a corotation resonance where |Aw| 
becomes very small and at a Lindblad resonance where 
V = 0. 

2.2. Driven Perturbations 

The influence of a small non-axisymmetric component 
of the stellar magnetic field can be studied by including in 
equation (lb) the small force due to this non-axisymmetry. 
At a given radial distance, the total vertical magnetic field 
B = B v + B Z consists of the vacuum component B v due to 
current flow inside the star and the induced field B % due 
to the current flow in the plasma outside the star. The 
vacuum field field has the general form 

B v =B v0 (r) + AB v {r,<j>,t) , 

AB" = B vl (r) exp(i<p—iCl^t)+B v2 (r) exp(2i</>— iQ*t) + ... , 

(21) 

where B v0 is the axisymmetric component, B vl <C B v0 
is the quadrupole component, and B v2 <C B v0 is the oc- 
tupole component. The quadrupole term can represent 
a non-centered but aligned dipole field in the star. The 
total magnetic force or acceleration is F = —V{[(B) + 
AB v ] 2 }/{8irp), (because (B-V)B=0), where (B) = (B 1 + 
B v ) and where the average is over (f>. Linearization gives 
SF = -V[(B}AB v (r,(f> 1 t)}/{ATrp). Thus wc have 

(SF r , SF^) = (SF r0 , SF^o) exp(im0 - tfi.t) , (22) 

where m = 1 or 2 and f2* is the angular rotation rate of the 
star. The calculation of §2.1 is modified slightly beginning 
with equations (7) which are here replaced by 

A6v r + BSv^ = f r + SF r , 

CSv^ + D5v r = U + SF^ , (23) 

where A, ..,D are the same as defined below equation (7). 

The forcing term SF gives rise to an additional contri- 
bution to pSv not included in §2.1. This contribution is 
1 
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where V = V/ (rp) and V is the Lindblad factor defined in 
equation (10b). Including this contribution we find that 
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response of the flow ip = */ (P) ' " is given by 

r(V) 1/2 



if" - U{r)<p 



Au 



V H, 



(25) 



5 



where U is defined in equation (17) and H is a known func- 
tion determined by AB" . We are interested in only the in- 
homogeneous solutions to this equation where U, u>, Aw, 
and T> are real. 

2.3. Driven- Modulated Perturbations 

Here we consider the case where the abovementioncd 
driving force SF is modulated by the flow perturbation 
^> (LR07). This modulation comes about naturally by in- 
cluding the magnetic field perturbation SB in the magnetic 
force F. That is, F = -V{[(B) + SB + AB V } 2 } / (87175). 
Linearization of this force gives a contribution F" 1 = 
-V[8BAB v ]/(4irp). The equations of §2.1 can be used 
to derive an exact formula relating SB to 'J', but here we 
use a simplified formula for one of the dominant terms, 
SB = B^iVL?^)" 1 , in order to limit the complexity of 
the equations. An equation for the driven modulated per- 
turbations is then obtained by the replacement in equa- 
tion (23) of SF -> SF m = ^*SF/(VLl), where SF is still 
given by equation (22) and is a known function determined 
by AB V . The use of \&* rather than \f is required due 
to our assumed dependence of <5F in equation (22). For 
the driven-modulated perturbations we obtain in place of 
equation (25) 



tp" — U(r)ip = i- 
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where we have used the fact that \& = (VY^tp. The hats 
over different quantities indicates that they are now oper- 
ators with ijj — > i[d/dt) and — > —i{d/d(f). The matrix 
M allows us to write equation (24) as H = M • SF. 

Equation (26) is a linear equation for tp so that in gen- 
eral 

if = ipo exp(— iu>t)+cpi exp(i<fi—iuit)+(p2 exp(2i<j)—iu)t)+.. , 

(27) 

where uj is real but undetermined at this point. Clearly, 
the time and angle dependence on both sides of equation 
(26) must match. 

For the case where there is a small quadrupole field com- 
ponent AB = B vl cxp(icf) — iu>t), equation (26) implies 

<f%-V , u (r)<po=[Op(<pt)] 0u (m = 0), 



Pi-tfi,u,(r)¥>i=[Opfo>5)] lu (m=l), 



(28) 



where Op stands for the linear operator on ip* on the right- 
hand side of equation (26) and the new subscripts indicate 
the m value and the lu value. Here, we necessarily have 

as found earlier by LR07. 

For the case where there is a small octupole field com- 
ponent AB — B v2 exp(2icf> — icvt), equation (26) implies 

V>'l-Ux,M^i= [Op(vi)] 1)W , (30) 

where we again have oj — 0*/2. 

In contrast with the case of driven perturbations of §2.2, 
equation (26) is linear in tp so that its magnitude is in- 
determinate. Further study is needed to determine the 
magnitude of tp in this case. One possibility is to gen- 
eralize equation (26) to include on the right-hand-side of 
the equation the inhomogeneous driving force oc SF as 



well the nonlinear force SF nl = -X7(\5B\ 2 )/{8irp). Ow- 
ing to the nonlinearity the driving force SF at frequency 
fi* can generate the 1/2 subharmonic oscillations of tp at 
the frequency f2»/2 described by equation (26). Analogous 
subharmonic generation is known in similar systems (e.g., 
Minorsky 1974). 

3. MODEL OF EQUILIBRIUM 

Figure 2 shows sample measured profiles of the midplane 
axial magnetic field (B), the azimuthal frequency of the 
disk matter =v$/ (27rr)] and its density (p) from MHD 
simulations by Romanova, Kulkarni, & Lovelace (2008) 
for globally stable case. The simulations were three di- 
mensional for an approximately axisymmetric case. These 
profiles motivate our choice of analytic functions to repre- 
sent these quantities. 

The fact that v^{r) decreases as r decreases close to 
the star is due to magnetic braking. A small twist of the 
star's magnetic field transports angular momentum of the 
disk matter to the star. That is, for z > there is a ver- 
tical flux of angular momentum — rB^B z /(AT:) > with 
\B§\ -C \B Z \ which transports the disk angular momen- 
tum to the star along the star's field lines (see Figure 1). 
This loss of angular momentum implies a mass accretion 
rate of the disk (outside of the region of the funnel flow) 
of Md — ~r 2 (B c f > )f l B z /(d£/dr) (= const) in the absence 
of viscosity (see Lovelace, Romanova, & Newman 1994), 
where the h— subscript indicates evaluation at the top sur- 
face of the disk, t — ru$ is the specific angular momentum, 
and d£/dr is positive for the considered profiles. The tur- 
bulent viscosity due to the magnetorotational instability 
(MRI) is absent in the region of the disk where ca > c s 
and/or dVl^/dr > (Balbus & Hawley 1998). Note how- 
ever that the disk equilibria discussed below neglect the 
accretion and have Ba, = 0. 




Fig. 2. — Midplane radial profiles of the magnetic field, the az- 
imuthal frequency [vj, = v l j > /(2nr)], the Keplerian frequency {vk)i 
and the density p from the MHD simulations. The vertical scale is 
for the azimuthal frequency of the disk matter; the star's frequency 
is = 300 Hz. The scalings of B z and p is arbitrary. 

We assume a pseudo-Newtonian gravitational potential 
§g = -GM«/(r — rs), where M* is the star's mass and 
r s = 2GM*/c 2 = AAA x 10 5 cm for a 1.4M Q star. The 
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angular velocity of a single paraticle is fix = 2-kvk = 
{GM*/[r(r - r s ) 2 }} 1/2 for r > 3r s . Near the star, the ro- 
tation frequency of the matter is modeled following LR07 

as 

v l + /(r) l + 

which is a first approximation as explained below. Here, z/* 
is the rotation frequency of the star and f(r) = exp[— (r — 
ro)/A] with ro the standoff distance of the boundary layer 
and A is its thickness. The azimuthal velocity of the mat- 
ter is U0 = 27rrz/J. Both ro and A are expected to depend 
on the accretion rate and the star's magnetic field. 
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Fig. 3. — Radial profiles of the magnetic field (B z ), the azimuthal 
frequency [kj, = v& / (2irr) . the Keplerian frequency (vx), and the 
density p of the model discussed in §3. The horizontal dotted line 
indicates the real part of the frequency v r of the unstable m = 1 
corotation mode discussed in §4.1. The vertical scale is for the az- 
imuthal frequency of the disk matter; the star's frequency is taken 
to be u f = 300 Hz. The plot scales for B z and p are arbitrary. For 
this case, rg/rg = 4 and A/rg = 0.1 in equation (25), e = 0.05 and 
po = 0.013 g cm -3 in equation (26), 7 = 5/3, and at r/rs = 6, 
c A /c s = 0.156. For this plot B z (r = ro) = 9.32 X 10 8 G. 



The radial force equilibrium for the axisymmetric flow 



is 



d 
dr 



P- 



52 
8tt 



(32) 



where p = pc 2 /j is the pressure in the disk and c s is the 
sound speed. The density profile is modeled as 



Po 



(1-e) exp(-0.05r/r s ) 

r+7w 



(33) 



with /(r) the same function as in equation (31) and e is a 
positive quantity much less than unity. The sound speed 
is modeled by choosing say c s k — O.lrSlx(r) and then 
letting 

c S K(r )f(r) c sK {r) 
1 + f(r) 1 + f(r) 

At large distances r 3> tq, c s jv$, = 0.1 corresponds to 
a disk half-thickness h w O.lr. Inside of ro the sound 
speed is assumed to be a constant. Using these expres- 
sions we first solve equation (32) for B 2 using equation 
(31) and neglecting the pressure gradient. We then go 



° + Sv^. The correction is small, 
<C 1. The radial epicyclic frequency 



back and obtain the correction <5z/^ to the azimuthal fre- 
quency needed to account for the pressure gradient; that 
is, 2pr(2-K) 2 v^ > 8v t j } — dp/dr so that the actual rotation 

frequency is v$ = v 
\6v^\ ~ (c s /« ) 2 
is calculated from Q r = 2-nv r with v 2 = r^ 3 'd(r A vi) / 'dr . 

Representative curves are shown in Figures 3 and 4. The 
value of the magnetic field at r/rs — 6 is arbitrary, but 
here it is chosen so that ca < c s which allows the MRI to 
grow which in turn gives rise to a turbulent viscosity in the 
disk (Balbus & Hawley 1998). The maximum of the disk 
frequency v$ is about 1149(4rs/ro) 1 ' 68 Hz for A/rg = 0.1, 
and it occurs at a distance about 0.26rs larger than ro- For 
A/rs — 0.2, the maximum of v$ is somewhat smaller and 
it occurs at a distance about 0.4rg larger than ro. 
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Fig. 4. — Radial dependences of the Alfven speed, Cj\, sound 
speed c s , and azimuthal velocity for the same conditions as in 
Figure 3. In this plot the speeds are in units of 10 10 cm/s. 

4. RESULTS 

We solve the equations of §2 using the equilibrium pro- 
files of §3. For this we use Maple V. 12 where we de- 
fine about 40 different functions of R = r/rs, for exam- 
ple, p(R),is^(R),v r (R),Auj(R),B(R),..,U(R). Some of 
the functions, for example, U , are complex for complex 
lo. The axisymmetric (m = 0) perturbation are found to 
be stable if the width of the boundary layer is not too 
narrow; that is, there is stability for A/rs ^ 0.02. In 
the following we consider the possible instability of modes 
with m — 1,2, .... 

4.1. Free Perturbations 

Figure 5 shows the radial dependence of the real part 
of the key function <7(r)^ r (r) for the same conditions as 
Figure 3. This function has a maximum at vr — 3.92rs 
indicated by the vertical arrow marked corotation reso- 
nance. The function changes sign at r^ = 4.07rs which 
is a Lindblad resonance where T (and T>) change signs. 
The maximum of 5R[(;(r)J r (r)] appears to be a general fea- 
ture for profiles similar to those in Figure 2. The depen- 
dence can be traced to the dependence of T>{r) (equa- 
tion 10b). For r decreasing significantly below ro, the 



7 



term —{dp* / dr) / (pL*) in T> becomes increasingly nega- 
tive. This is due to the radial dependence of the magnetic 
pressure and the small density. Note is negative in this 
region. On the other hand for r increasing from ~ ro, the 
positive contribution of fl 2 , begins to dominate. The com- 
bination of these dependences gives a $l[g(r)F(r)] profile 
with a maximum at a distance inside the Lindblad reso- 
nance as shown in Figure 5. 



1 
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Fig. 5. — Radial dependence of the real part of the key function 
gT for for the same conditions as for Figure 3. The imaginary part 
of the function is small compared with the real part. The function 
has only a weak dependence on w which for this plot is w ~ 850 + i89 
Hz. We find maximum instability for SR[Ao;(r)] = at the radial 
location of the maximum of ^(gj 7 ) indicated by the vertical arrow. 
The other arrow shows a Lindblad resonance where $t(F) changes 
sign. The vertical scale is arbitrary. 

Figure 6 shows the real part of the effective potential 
U(r) for the same case as Figure 5. The real part of the 
frequency is chosen to give SR[Aw(r)] = at rji/rg = 3.92 
of the maximum of $l[g{r)!F(r)\ because this give the max- 
imum growth rate. For m = 1 this gives v r = uj r j2ix = 850 
Hz. Clearly this frequency must be larger than rotation 
frequency of the star . 

The depth of the potential well shown in Figure 6 in- 
creases as the imaginary part of the frequency uji = 5 (a;) > 
(the growth rate) decreases. We use a WKBJ treatment 
with <p cx fc- 1 / 2 exp(±i f drk) and k = [-^(C/)] 1 / 2 . 
The imaginary part of U is small compared with the real 
part. The allowed values of Wj are then calculated using 
the Bohr-Sommerfeld quantization condition, 



dr k 



0,1,2,. 



(35) 



where k = — 9?(f7), and n n , & r out are the radii where 
R(Z7) = 0. 

For the case of Figure 6, r in /rs — 3.83 and r out /rs — 
3.97. The largest growth rate u>i corresponds to n = 
and for the case of Figure 6 with m = 1 this gives 
Vi = u>i /2tt = 89 Hz so that u>i/u) r as 10%. For modes 
with m > 2 also grow with similar cjj/w r values but as 
discussed below these modes do not give time variations 
in the total flux from the source. 

For a more gradual boundary layer with A = 0.2 and 
m = 1, we find Sft[Aw(r)] = for v r = 750 Hz at the ra- 



dius r/rs — 3.83 of the maximum of ^R[g(r)J 7 (r)]. From 
equation (35) we find v. L = 120 Hz. For A increasing from 
0.2, we find that v r and rn continue to decrease gradually 
and vs also decreases. 
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Fig. 6. — Radial dependence of the real part of the effective poten- 
tial 3t[f/(r)] (multiplied by r|) for the same conditions as for Figure 
5. The imaginary part of U is much smaller than the real part. 

Owing to the perturbation, the surface temperature of 
the disk is 

T(r, <f>, t) = T + 5i[(5Ti exp(i</> - iu x t) 

+ ST 2 (r) exp(2i0 - iuj 2 t) + . .] , (36) 

where Tq(t) is the unperturbed temperature, STi^(r) -C 
To are the amplitudes of the m = 1, 2 corotation modes, 
and u>i.2 are their frequencies. We neglect the difference 
between the midplane and the surface temperature of the 
disk. The corresponding flux density from the disk sur- 
face is proportional to S(r, cj>, t) - T 4 + K[4T 3 (5Ti exp(i0- 
iuit) + 4T$6T 2 exp(2i(f> - ioj 2 t) + ..] . The total flux for a 
face-on disk, L ~ J rdrd<pS, is independent of time: the 
<p— integration annihilates the terms dependent on (j>. For 
a more general disk orientation, with the disk angular mo- 
mentum tilted say towards the line of sight by an angle l, 
the Doppler effect due to the disk rotation gives a boost 
to the frequency for say <j> — and a decrement for <f> = it. 
This corresponds to multiplying S by the Doppler factor 
D(r,(j)) = [1 + e(r) cos(</>)] 4 w 1 + 4ecos(</>), where e = 
[v<t,(r)/c] sin(t) and e 2 1. For the case of Figure 6 where 
r = 3.92r s , v^/c — 0.289. Consequently, there is a contri- 
bution to the source flux SL ~ J rdrd(f>D[r,(f))S{r,(j),t) ~ 
167r5R[/ rdrT^eSTi(r) exp(— iwtt)] . We interpret this fre- 
quency as the upper frequency component of the twin 
QPOs as argued earlier by LR07. Note that the higher 
order terms ST 2 , ST3, ... give no contribution to the total 
flux due to the 4>— integration. 



4.2. Driven Perturbations 

Here we consider the driven perturbations of the flow 
which result from a quadrupole or octupole component 
of the star's magnetic field as discussed in §2.2. Inspec- 
tion of equation (25) reveals that the radial locations of 
the Lindblad resonances where T> (or T>) vanishes are very 
important for excitation of the flow (LR07). At such a 
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resonance the right-hand-side equation (25) - the driving 

term - is proportional to T> (2?)~ 3 / 2 . 

We first consider the quadrupole field case B vl where 
m = 1 for v r = cj r /2w = 300 Hz and io i = 0. We find 
that T> goes through zero at one radius, = 4.07rs for 
this case. Near this radius we find T> ss a(x — bx 2 ), where 
x = (r — r/,)/rj, and a, b > are constants. We are in- 
terested as explained below to determine the amplitude of 
this driven mode near the vicinity of the corotation mode 
at rfl = 3.92rs. For this purpose we develop an approxi- 
mate solution to equation (25) for x 2 <C 1 by retaining only 
the most singular term in the effective potential which is 
U « (3/4) (v'/V) 2 . The right-hand side of the equation 

can be approximated as KT> (2?)~ 3 / 2 , where K oc B vl . 
Thus equation (25) simplifies to 

„a.__ K(l 

dx 2 



K 2 ip 



2bx) 



(37) 



(x -bx 2 ) 3 / 2 ' 

where n 2 = (3/4)[(l - 2bx)/(x - bx 2 )} 2 . The inhomoge- 
neous solution to this equation for x 2 <C f is 

Although ip diverges at x = 0, note that 4" = Sp*/p is a 
constant. For —x increasing we find that the dependence 



of n changes from 



(3/4)a;- 2 to k 2 



38. The 



small value of Kg means that driven mode has a signifi- 
cant amplitude at the distance r# where the free mode is 
excited. 

With both the free perturbation and the quadrupole 
driven perturbation present we have 

T(r, <p, i) = T +$l[6T{(r) exp(^ - iu) X t) 

+ 5T?(r)exp(i^-iQ„t)] . (39) 

The associated flux density S ~ T 4 is 



7^4 
-f n 



ifLi) 



- „ -r 4T 3 5R [ST{ exp(i(j> - iuj x t) + ST? 

+6T 2 ^R{ST{ ST?* expH(o;i - ft*)i]} + (40) 

where the ellipsis denotes terms of the form 0[exp(±2i(f>)] 
and 0[exp(±4i<^))] which do not cause variations in the to- 
tal flux. It follows from this equation that there are three 
QPO components for a generally oriented disk. Two of the 
components arise from the above-mentioned Doppler boost 
acting first on the term ST( exp(i^) — iw\t), which gives a 
frequency ui\ component in the total flux, and second on 
the term ST? exp(icf> — zft*i), which gives a frequency ft, 
(the star's rotation frequency) in the total flux. The third 
component arises for the final term in equation (37) and 
has a frequency u>i — ft*. The component at ft* is how- 
ever usually absent (van der Klis 2006) so that we do not 
consider this case further. 

For the case of an octupole field component B v2 where 
m = 2 we find two Lindblad resonances for o; r /27r = 
ft*/27r = 300 Hz and oji — as shown in Figure 7. One res- 
onance is at r^i/rs = 4.11 and the other at rLo/fs = 4.28. 
The driven motion at ru, is important here because this 
is close to the radius of the unstable free perturbation rfj. 
For this case 

T(r, t) = T + &[5T{(r) exp(# - iuj x t) 

+ 5T?{r) exp(2i0 - ift*t)] . (41) 



The associated flux density is 
T 4 + 4T 3 3?[(5T/ exp{icj) - iwit) + ST? exp(2^ - ift**)] 
+6T 2s R{5T{ ST?* exp[-i(j) - - ft*)i]} + ... . (42) 

In this case we have just two frequency components: The 
first is at the frequency u>x of the unstable free per- 
turbation as a result of the Doppler boost acting on 
the term 8T[ ex-p(icf> — iu>xb). The second is at the fre- 
quency uji — ft* as a result of the Doppler boost act- 
ing on the term STiST?* exp[-i(f> - i{w x - ft*)t]. The 
Doppler boost acting on the remaining term in equation 
(39), ST!? exp(2i(f) — itt*t), causes no variation in the total 
flux. 
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Fig. 7. — Radial dependence of the Lindblad factor T> (arbi- 
trary scale) and the effective potential U (multiplied by r|) for 
uv/27T = 300 Hz, u>i = 0, and m = 2 and other conditions the same 
as in Figure 5. 

4.3. Driven Modulated Perturbations 

Here we consider the driven modulated perturbations of 
the flow which result from the octupole component of the 
star's magnetic field B v2 discussed in §2.3. Two Lindblad 
resonances are again found at radii similar to the case of 
Figure 7. The inner radius rn is important here because it 
is close to the radius rji of the unstable free perturbation. 
For the driven modulated octupole case we have 

T(r, cf>, t) = To + R{ST((r) exp(i<f> - iuit) 

+ ST? m (r)exp[2i<j>-i(n*/2)t}} . (43) 

The associated flux density is 

T 4 + 4T 3 3?{ 5T( exp{i<j) - iuit) + ST? m exp[2i^ - i (ft* /2)t)} 
+6Tffi{5T{ST? m * exp[-i<p-i(ui - ft*/2)t]} + ... . (44) 

In this case we again have just two frequency compo- 
nents: The first is at the frequency u>i of the unstable 
free perturbation as a result of the Doppler boost act- 
ing on the term ST[ exp(i(j> — iuiit). The second is at the 
frequency wi - fl,/2 as a result of the Doppler boost act- 
ing on the term ST(5T? m * exp[-i<f) - i(u) X - ft*/2)i]. The 
Doppler boost acting on the remaining term in equation 
(39), ST? m exp[2i(j> — i(ft*/2)t], causes no variation in the 
total flux. As mentioned in §2.3 the present theory does 
not predict the value of 5T? m . 
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5. NONLINEAR EFFECT OF UNSTABLE WAVE 

The considered equilibrium disk (§3) does not include 
accretion. However, the unstable Rossby mode discussed 
in §4.1 gives rise to accretion in the vicinity of t\r, 

M Ro = -2nr f dz $t{5p5u*) , (45a) 

J-h 

where h is the half-thickness of the disk. We evaluate 
this by taking the most singular terms in Auj in the limit 
where this quantity almost vanishes at r — tr. This gives 
dp = 2k (j> J r ^/(Au>L^), and 6u r = -2ik^T^ j p. Thus, 

8wj/ifc5|.F| 2 |*| 2 

MR ° = IA 121- < ( 456 ) 

where u>i > is the growth rate and < 0. Thus, the 
instability acts increase the mass accretion rate, Mr q > 0. 

The unstable Rossby mode also gives an inflow of angu- 
lar momentum, 

r-h 



-2-KT 



dz 



Using the fact that pSu^, - 
find 

F Ro = M Ko lQ , 



$l(5p5u*)ru<i, + p$t{5u*r5u<{,) 

(46a) 

-k^^p^/ipfl^AwLt), we 



(466) 



where Q = (l + ?i| f /u|)/2 and uk = tQ,k is the Keplerian 
velocity. For the case of Figure 6, Q = 1.884. Because 
Q > 1, the Rossby mode transports inward more specific 
angular momentum I than exists in the equilibrium. Con- 
sequently, a nonlinear effect of the Rossby mode is to make 
the positive slope of £(r) smaller in the vicinity of tr. A 
sufficiently large reduction of dljdr may cause the growth 
rate to decrease. 

6. CONCLUSIONS 

We have investigated three modes of accretion disks 
around rotating magnetized neutron stars which may ex- 
plain the separations of the kilo-Hertz quasi-periodic oscil- 
lations (QPO) seen in low mass X-ray binaries. This work 
is a continuation of the earlier work by LR07 where these 
modes were identified. We develop the theory of compress- 
ible, non-barotropic MHD perturbations of an axisymmet- 
ric equilbrium flow with d/dz = 0. We assume that there 
is a maximum in the angular velocity Jl</,(r) of the accret- 
ing material larger than the angular velocity of the star 
£1*, and that the fluid is in approximately circular motion 
near this maximum rather than moving rapidly towards 
the star or out of the disk plane into funnel flows. MHD 
simulations by Romanova et al. (2002, 2008), Long et al. 
(2005), and Kulkarni & Romanova (2008) show this type 
of flow and f2^(r) profile. 

The first mode we find is a Rossby wave instability or 
unstable corotation mode which is radially trapped in the 
vicinity of the maximum of a key function g(r)T{r) at 
tr. We derive a Schrodinger type equation for the per- 
turbation, <p" = U(r)(p, where U(r) is the effective poten- 
tial. This instability is analogous to that found earlier by 
Lovelace & Hohlfeld (1978), Lovelace et al. (1999), and 



Li et al. (2000) and in simulations by Li et al. (2001). 
The real part of the angular frequency of the mode is 
iv r = mVt^rR), where m — 1,2... is the azimuthal mode 
number. The imaginary part of the frequency u>i (the 
growth rate) is determined by a Bohr-Sommcrfcld quan- 
tization of the perturbation in the effective potential. We 
argue that for generally oriented disk, the Doppler boost of 
the disk surface emission from the unstable m — 1 mode 
will give periodic variations in the total flux with angu- 
lar frequency uj r which we suggest is the higher frequency 
component of the twin QPOs as proposed by LR07 The 
saturation of the growth of this mode remains to be ana- 
lyzed in detail in future work. 

The second mode, is a mode driven by the rotating, 
non-axisymmetric quadrupolc [~ cxp(z0)] and/or octupolc 
[~ exp(2z(/>)] components of the star's magnetic field. It 
has an angular frequency equal to the star's angular ro- 
tation rate £7*. This mode is strongly excited near the 
radius of the Lindblad resonance which is slightly outside 
the radius of the maximum of g T . When both the first and 
second modes are present the nonlinearity of the emission 
will in general give a product term with angular frequency 
uj r — £1*. For a quadrupole field taking into account the 
Doppler boost, we find that the total flux has periodic 
variations at three frequencies, u) r , £1*, and uo r — £7*. How- 
ever, the frequency £1* is usually not observed (van der 
Klis 2006). The situation is different for an octupole field 
component again taking into account the Doppler boost. 
In this case the total flux has periodic variations at two 
frequencies, u r the upper QPO and ui r — £1* the lower 
QPO. 

The third mode arises from the interaction of the flow 
perturbation with the rotating non-axisymmetric compo- 
nents of the star's magnetic field. This interaction arises 
naturally owing to the fact that the magnetic force is 
proportional to the gradient of the square of the mag- 
netic field. We derive a linear differential equation for 
this driven modulated perturbation. We find that the an- 
gular frequency of the perturbation is £l*/2 both for the 
quadrupole and the octupole field components. In the case 
of the octupole field component and taking into account 
the Doppler boost, we find that the total flux has peri- 
odic variations at two frequencies, iv r the upper QPO and 
Lo r — Q*/2 the lower QPO. Because the equation for the 
driven modulated perturbations is linear, the amplitude 
of the motion is indeterminate. One possibility is that the 
0*/2 motion is excited nonlincarly by the rotating non- 
axisymmetric field which has frequency fi*. This will be 
discussed in a future work. Thus present theory does not 
determine whether the lower QPO frequency is u> r — £1* or 
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